function [gamma] = gammln(x)
% gammln -
% from Press et al., p. 214 
% used for t-test etc.
% coded P. Manis, 12/2/98
%
cof = [76.18009172947146, -86.50532032941677, 24.01409824083091, ...
      -1.231739572450155, 0.1208650973866179e-2, -0.5395239384953e-5];
y = x;
tmp = x+5.5;
tmp = tmp - (x+0.5)*log(tmp);
ser = 1.000000000190015;
for j=1:6
   y=y+1;
   ser = ser + cof(j)/y;
end
gamma = -tmp + log(2.5066282746310005*ser/x);
return

   